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A Gamma-ray Pulsar Timing Array Constrains the 
Nanohertz Gravitational Wave Background 
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After large galaxies merge, their central supermassive black holes are expected 
to form binary systems whose orbital motion generates a gravitational wave 
background (GWB) at nanohertz frequencies. Searches for this background 
utilize pulsar timing arrays, which perform long-term monitoring of millisec- 
ond pulsars (MSPs) at radio wavelengths. We use 12.5 years of Fermi Large 
Area Telescope data to form a gamma-ray pulsar timing array. Results from 
35 bright gamma-ray pulsars place a 95% credible limit on the GWB charac- 


teristic strain of 1.0 x 10714 at 1 yr ^!, which scales as the observing time span 


—13/6 
obs 


t . This direct measurement provides an independent probe of the GWB 


while offering a check on radio noise models. 
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Main Text 


Pulsars are spinning neutron stars that emit beams of broadband radiation from radio to gamma- 
ray wavelengths that appear to pulse as they periodically sweep across the line of sight to 
Earth (/). Millisecond pulsars (MSPs) spin at hundreds of hertz and pulse with sufficient reg- 
ularity to function as celestial clocks distributed across the sky and throughout the Galaxy. 
Timing of individual MSPs using radio telescopes has been used to test General Relativity and 
alternate theories of gravity (2). Long-term monitoring campaigns of ensembles of MSPs are 
used to search for low-frequency gravitational waves, expected from supermassive black hole 
(SMBH) binaries that are predicted to exist at the centers of galaxies that have undergone merg- 
ers. General Relativity predicts that a circular binary with orbital frequency f /2 will emit GW 
with frequency f and amplitude œ f?/? (3), and when SMBH binaries have an orbital separation 
of ~0.01 parsecs (pc, ~2000 astronomical units), the orbits decay primarily through GW emis- 
sion. Because of this link between GW frequency and amplitude, the superposition of GWs 
from many SMBH binaries throughout the Universe is predicted to build up a GW background 


(GWB) with a characteristic GW strain ^, following a power law in frequency (4), 


hf) dios (4) (1) 


yr 
The spectral index a is predicted to be —2/3 for GW-driven binary inspirals, while the dimen- 
sionless strain amplitude A, incorporates the growth, masses, and merger rates of SMBHs. If 
SMBHs do not rapidly migrate to the centers of newly-merged galaxies, there will be relatively 
fewer wide binaries, reducing the GW power at low frequencies. Thus the measured GWB 
is expected to carry information about the distribution of SMBH masses and the dynamical 
evolution of SMBH binary systems (5). 


This GWB can be detected with ensembles of MSPs—known as pulsar timing arrays (PTAs) 


(6, 7)—by monitoring the times of arrival (TOAs) of the steady pulses from each pulsar, which 
arrive earlier or later than expected due to the spacetime perturbations. Because the GWB is 
expected to be a sum of many individual sources, the induced TOA variations are random and 
differ for each pulsar, but have a common spectrum of power spectral densities, P( f): 


A? -r 
PE (4) yr, (2) 


with spectral index I = 3—2a = 13/3 for SMBHs (4). This functional form has more power at 
low frequencies so is referred to as a red spectrum. For observations taken at an approximately 
fixed location (Earth), the GWB is expected to produce a signature quadrupolar pattern of TOA 
variations, known as the Hellings-Downs correlation (8). 

Because the expected quadrupolar correlations are only about 10% of the total signal, the 
GWB is predicted to first appear as a set of independent signals from each pulsar whose power 
spectra are all consistent with Equation 2, with the quadrupolar distribution only becoming 
evident in more sensitive observations. Radio PTAs have reported a red spectrum process with 
modest statistical significance (9-/2), but no Hellings-Downs correlation has been found. These 


results could be compatible with a = —2/3 and Ag, ~ 2 — 3 x 10 P at 1yr^! 


; see Figure 
1. This would be consistent with some predictions for the GWB (5), but because no spatial 
correlations have been detected, it could have other origins. 

A potential alternative explanation for this signal is spin noise, approximately power-law 
red noise intrinsic to each pulsar, with some MSPs observed to have a spin noise spectral index 
(1°) of 2-7 (13, 14). Possible physical origins for spin noise include turbulence in the neutron 
star interior (/5) and systematic variations in the magnetic field and co-rotating plasma which 
govern the rotational energy loss of the pulsar (76). Pulsars that have spin noise spectra with 


similar shapes but different amplitudes—inconsistent with a GWB—could masquerade as a 


common mode signal without a Hellings-Downs correlation (/0). 
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Figure 1: Constraints on the gravitational wave background from radio and gamma-ray 
PTAs. The inferred constraints on the GWB amplitude at 1 yr ! (Agy», data sources in Table 
S7) are plotted as a function of publication date and assume a = —2/3, as predicted for the 
superposition of gravitational waves from merging black holes. Colored symbols, indicated in 
the legend, correspond to each of the PTAs. Upper limits at 95% confidence are shown with 
arrows while amplitude ranges indicate detections of a common noise process. The Fermi-LAT 
95% upper limit, 1.0 x 10714, uses data obtained through January, 2021 and is shown as a red 
star at an approximate publication date April, 2022. The dashed red line indicates the expected 


scaling as the limit as a function of time. 


Another potential noise source for radio pulsar timing arrays is the frequency-dependent 
effect of radio propagation through plasma, including the solar wind and the ionized interstellar 


medium (IISM). Pulsed radio emission at frequency v is delayed by time Tpm 


415 © DM cm? si ( V 9 (3) 
TpM = 4.15m ———— 
EM j pe GHz/ ^' 


where DM is the dispersion measure, equal to the total electron column density. The DM of a 
pulsar can vary with time, due to the relative motions of Earth and the pulsar. Correcting for this 
effect requires repeated measurements using multi-frequency radio observations and the intro- 
duction of many additional degrees of freedom to timing models. Because the propagation paths 
of radio waves through the IISM depend on v, the DM itself is frequency-dependent (/7, 16), 
so some of this delay is intrinsically unmeasurable. Other propagation effects include a broad- 
ening of the pulse that can only be corrected for bright pulsars, with some components also 
being unmeasurable (/9). Because the IISM is turbulent, these uncorrected delays introduce 
additional red noise to radio pulsar timing data. The variable solar wind introduces similar 
dispersive delays which can in principle be measured like DM variations but are only partially 
included in current models (20). Due to the wide angular extent of the solar wind, uncorrected 
delays would be correlated amongst pulsars. As with spin noise, IISM-induced noise with sim- 
ilar spectra could mimic a GWB signal. Predicted noise amplitudes are similar to the expected 
GWB signal, but these predictions rely on assumptions about the turbulent spectra of the IISM 
which are poorly constrained by data (19). See (21) for further discussion of the modeling and 
impact of noise. 

Gamma-ray observations offer a potentially complementary approach: the much higher pho- 
ton frequency means that the effects of the IISM and solar wind are negligible. The Large Area 
Telescope (LAT) (22), on the Fermi Gamma-ray Space Telescope, is sensitive to GeV gamma- 


ray photons emitted by MSPs. Its 2.4 steradian field-of-view performs a continuous survey, 


covering the full sky every 2 orbits (~3 hr). Its GPS clock records photon arrival times with 
<300 ns precision (23), enabling pulsar timing. Analyses of the LAT survey data have de- 
tected 127 (21, 24) of the over 400 known MSPs in the Milky Way. The large MSP sample, 
long observing span, and instrumental stability enable a gamma-ray pulsar timing array whose 
characterization of spin noise and a potential GWB signal is free from IISM effects. 

Using the 35 brightest and most stable y-ray MSPs and 12.5 yr of Fermi-LAT data, we 
searched for the GWB using two different techniques (2/). First, we implemented a coherent 
photon-by-photon analysis which retains <1 us resolution. Second, for analysis with the estab- 
lished software used for radio data analysis, we directly measured TOAs from the LAT data (25). 
Because the TOA estimation procedure requires averaging up to one year of data, this method 
loses sensitivity to shorter-timescale signals, and only 29 of the 35 pulsars are suitable. 

For each pulsar, we searched for spin noise and derived an upper limit on A4, using the 
photon-by-photon method and with two TOA-based software packages, TEMPONEST (26) and 
ENTERPRISE (27). None of the pulsars show evidence for spin noise (27) and the three different 
methods provide consistent results for each pulsar (Figure 2) except for three cases (21). 

Three of the pulsars in our sample have spin noise measurements from radio PTAs. Using 
the power spectral indices I' measured from the radio timing data, we calculated 95% upper 
limits on spin noise amplitudes from the gamma-ray data. Our limits are below the previously 
measured values for PSR J0030+0451 (10% of the measured value) and PSR J19394-2134 (60— 
77096) but are unconstraining for PSR J0613—0200. This discrepancy might indicate contami- 
nation by residual IISM effects on the radio-based spin noise and GWB signal measurements. 

We combined the single pulsars into a pulsar timing array and estimated A, limits under 
a variety of scenarios, including marginalization over possible spin noise and uncertainties in 
the position of Earth relative to the Solar System barycenter, and both excluding and including 


the quadrupolar spatial correlations expected under General Relativity (27). The resulting rep- 
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Figure 2: Comparison between A wp measurements from each pulsar using three analysis 
methods. Data points indicate the limits on an a = —2/3 GWB for 35 MSPs computed with 
three methods: two TOA-based codes TEMPONEST (orange stars) and ENTERPRISE (blue cir- 
cles) are shown as a function of the limit from a photon-by-photon analysis (x axis). The dashed 
line indicates equality between the results of the TOA-based and photon-by-photon methods. 
Six pulsars (purple triangles) have only a photon-based analysis so are plotted arbitrarily at zero 


on the y axis. The three labeled pulsars are outliers (27). 


resentative 95% confidence limit is Agw < 1.0 x 10- (Figure 1), a factor of 3-5 higher than 
the red spectrum process detected by radio PTAs. 

For an idealized PTA, when a potential GWB signal is weak compared to other noise, the 
signal-to-noise ratio grows proportionally to AS x th (28, 29), with tops the observing time 


obs 


span and I = 13/3 for SMBHs as in Equation 2. This means that upper limits on Ayyw improve 


—13/6 


following the relation Agyy o top. 


. On the other hand, if the signal that terrestrial PTAs 
are currently detecting does arise from the GWB, then these PTAs are now in the strong signal 
regime and their sensitivity will improve slowly (cc i. ^). The differing time scalings and noise 
sources allow the gamma-ray PTA data to distinguish residual IISM variations from a potential 
GWB signal. 

The Fermi PTA data have an essentially constant experimental setup: the data are almost 
uninterrupted and calibrations have been constant for the full 12.5 year dataset. Gamma-ray 
data are potentially less subject to astrophysical effects such as changes in the radio pulse shape 
(21). This stability is particularly useful for probing GWs with frequencies below 0.1 yr !. 
Such low frequencies are predicted to constrain the spectral shape of the GWB which contains 
information about the physical sources (5). 

There are other potential sources of power-law GWBs with different spectral indices, a, 
such as a = —1 for relic GWs originating during scale-invariant inflation in the early Universe 
(30). Decay of (hypothetical) cosmic strings could also produce power-law spectra under a 
variety of scenarios (3/). To constrain such sources, we computed corresponding 95% upper 
limits on A,» at different values of o (Figure 3). Yet other models are not well described by 
power laws, but predict the largest signal in or near the PTA band (32, 35). 

To summarize, we have used the Fermi LAT dataset to construct a gamma-ray PTA. This 


provides an independent method to search for signals detected by radio PTAs; unlike the radio 


PTAs, it is free from the effects of the IISM. Most of the pulsars are amenable to the TOA- 
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Figure 3: Gamma-ray constraints on different types of GWB sources. GWB amplitudes 
Aswb for assumed spectral indices a in the shaded region are excluded with 95% confidence. 
The symbols indicate the values of o expected for SMBH binaries (red star, our fiducial result), 
gravitational waves generated during cosmic inflation (green triangle) and from hypothetical 


cosmic strings (blue circle). 


based approach, and the resulting datasets are small compared to those of radio PTAs, enabling 


analysis alongside radio PTA data with little additional computational burden. 
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Materials and Methods 


Pulsar Timing using Radio and Gamma-ray Observations 


The spin phase at time t, (t), increases by 1 each time a pulsar rotates, and because the pul- 
sation mechanism is fixed to the star, ¢(t) can be inferred by observing pulses. Pulsar timing 
is the measurement of these pulse arrival times (TOAs) and comparison with a timing model 
that predicts @(t, A). The astrophysical parameters A leave a characteristic imprint on the tim- 
ing residuals between data and model: a spin frequency error dv produces linear residuals, dv 
leaves quadratic residuals, a position error induces an annual sinusoid, etc. (35). The goal of 
pulsar timing is to measure À and thus characterize a wide range of astrophysical phenomena 
and constrain fundamental physics (2, 36—36). 

During a typical —0.1-1 hr observation with a radio telescope, light collected from the 
antenna is transduced, amplified, digitized, dedispersed, and filtered into frequency channels. 
Pulses are stacked into a single pulse profile by folding the data at the instantaneous pulsar spin 
period. The observation time is recorded using a precise clock, and the offset of the recorded 
pulse to this time yields the pulse TOA. The uncertainty of a TOA can be estimated from the 
random noise in the pulse profile, which arises from electronics noise, background radiation 
from astrophysical sources, and terrestrial spillover into the antenna. The resulting uncertainty 
is Gaussian, so a typical radio timing analysis optimizes the timing model parameters and as- 
sesses goodness-of-fit by minimizing x?(A). 

The Fermi-LAT, on the other hand, collects individual gamma rays. The arrival time of the 
ith photon, t;, is recorded with ~300 ns precision, but aside from its energy, a photon carries no 
additional information: it could be from any pulsar phase ¢ or even from a background source, 
so these t; cannot be interpreted as TOAs. In some cases, histograms in @ can be built up over a 
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long enough time—hours for bright pulsars, years for the faintest—that a pulse profile emerges. 
Analogously to folded radio profiles, comparing these histograms to the assumed template / (4) 
can yield a TOA (39), and x?(A) methods may be applicable. (Radio and gamma-ray pulse 
profiles differ, but we use f (4) for the pulse profile in the relevant band.) 

For many of the faint MSPs used in this work, there is insufficient integration time to build 
up a gamma-ray pulse profile and estimate a TOA. For instance, resolving the annual sinu- 
soidal residuals from a position error requires sampling of at least 2 TOAs per year, and prefer- 
ably faster to avoid a systematic error. Instead, we can use a timing model to evaluate the 
phase @(t;, A) at each individual photon time and then gauge the agreement of the resulting 
distribution with an assumed template f(@) using the likelihood (Equation S1). Because the 
LAT has a broad, energy-dependent angular resolution, photons from different sources over- 
lap and we must also account for the background, which we do by computing the probability 
weight 0 < w; < 1 that the 7th photon originates from the pulsar (40-42). Using a normalized 
( T f(è) do = 1) pulse profile model, the Poisson likelihood for the data, Laata, is 


log Laata(A) = log (wi flt A) - a - w;)). (S1) 


By maximizing Laata(A), we obtain optimal estimates for parameters A while preserving the 
<1 us resolution of the LAT. 

Unlike in the x?(A) radio case, there are no residuals because there are no direct phase 
measurements (TOAs) with which to compare the model. This makes it more challenging to 
assess the goodness of fit and select the most favored models (43). Our GWB analysis (below) 
uses both TOA-based and photon-by-photon (likelihood) methods. We view the latter as more 
fundamental, but the former allows consistency checks with existing methods used for radio 
PTAs. 


Noise in Pulsar Timing Array Data 


PTAs (6, 7) extend pulsar timing to ensembles (arrays) of pulsars to search for correlated signals 
such as low-frequency gravitational waves. These subtle signatures can be obscured by noise, 
and many noise sources must be accounted for in radio PTA data to reach a typical amplitude 
in the timing residuals of «100ns. To place the Fermi PTA results in context, we give an 
overview of the radio PTA noise budget below, drawing distinctions to gamma-ray observations 
where appropriate. These main noise sources are summarized in Table S1, and in (44). 
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Table S1: Sources of noise in radio and gamma-ray pulsar timing array data. The list is 
incomplete and qualitative: for each noise source, we attempt to give an estimate of the impor- 
tance (amplitude in residuals, or difficulty of complete mitigation) and complexity (degrees of 
freedom (d.o.f.) in existing models for a typical pulsar). Dashes indicate that the entry is not 
applicable, and a question mark that its impact is unknown. 


Radio Gamma ray 
Noise Source Impact d.o.f. Impact | d.o.f. Note 
White Noise 
Measurement moderate — | major — | Sensitivity is major limiting factor for gamma rays. 
RFI minor ?|- — | RFI varies widely between observing systems. 
Calibration minor ?!— — | Affects certain pulsars/observing systems. 
Jitter moderate 10s | — — | Jitter affects high signal-to-noise observations. 
Red Noise 
DM variation major 100s | — — | DM(t) drives radio PTA observing strategies. 
Solar wind moderate ~10s | — — | Solar wind mitigation is poorly supported. 
Scattering moderate 100s | — — | Affects some pulsars/low radio frequencies. 
Pulse variability | moderate 0-10s | — ? | No gamma-ray MSP pulse profile changes known. 
Discontinuities | moderate 10s | — — | LAT data are continuous, not a general property. 
Spin noise major 10s-100s | major 10s | Fewer d.o.f. needed for less precise LAT data. 
White Noise 


White noise («1 day time scales) cannot produce GWB-like signals, but it must be modeled 
to estimate the precision of timing model parameters. Most analyses use the parameters EFAC 
and EQUAD (45), which modify the measurement uncertainty c; of a TOA according to o? + 
EFAC?o? -- EQUAD?. Additional ECORR parameters can be used to describe correlated noise, 
e.g. jitter. Because the noise sources depend on observing frequency, pulsar, and observing 


system, tens to hundreds of such parameters may be needed to represent the excess white noise 
in a radio PTA data set. The mapping of these parameters onto the actual noise sources may be 
imperfect, leaving un- or over-modeled white noise in the residuals. 

For any relevant PTA pulsar, gamma-ray measurement uncertainties are dramatically larger, 
and poor sensitivity is the major weakness of gamma-ray pulsar timing. 


Jitter Intrinsic pulse shape variations and the finite number of pulses received in an observa- 
tion cause a TOA bias. This source of noise is called jitter (46, 47). It requires a dedicated noise 
model because its impact depends on pulsar brightness, which varies due to IISM scintillation 
(below), and because it is correlated over the observing bandwidth (46). Because <1 photon 
is received from each pulse on average, jitter is fundamental to gamma-ray observations and is 
accounted for in Poisson statistics. 
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Radio-frequency interference (RFI) RFI affects all radio observatories but is less severe for 
those in remote areas. Strong RFI can render an observation unusable. Fainter RFI can often 
be characterized as narrowband—contaminating portions of the observing band, intermittently 
or continuously—or impulsive—affecting much of the band briefly. RFI biases TOA measure- 
ments (48): impulsive RFI alters the observed pulse profile while narrowband RFI changes the 
effective observing frequency, inducing residuals from DM(t) corrections (below). RFI noise 
is typically white, but a changing RFI environment can introduce biases on longer time scales. 
There is no analogous interference in Fermi-LAT data. 


Polarization Pulsars are highly polarized radio sources (49), so recording a stable pulse pro- 
file requires accurate polarization calibration. Since most radio telescopes use an altitude- 
azimuth mount, the measured pulsar polarization angle depends on the source elevation and 
must be corrected to obtain the intrinsic value. If a radio receiver has non-zero cross-polarization, 
this alters the pulse profile (50). Using bright pulsars to characterize the receiver response as a 
function of parallactic angle (5/) can eliminate much of the apparent TOA variation. As with 
RFI, observational bias can couple this white noise to longer time scales. Gamma-ray pulsar 
polarization is unknown, but Fermi-LAT is very insensitive to photon polarization (52). 


Red Noise 


Red noise operates on longer time scales and includes signals with power spectral densities 
similar to that of the GWB, so it must be modeled. 


Spin Noise Spin noise, also called timing noise or intrinsic noise, is well known in young, 
slowly spinning pulsars (53), but is also present in MSPs with an amplitude that is related to 
the pulsar spin-down power (54). This intrinsic noise is generally adequately modeled with a 
power-law spectrum of the same form as Equation 2, and represented in data either as long- 
term correlations between the TOAs, with a non-diagonal covariance matrix (55, 56), or using a 
Fourier expansion with the coefficients constrained to follow the assumed power spectrum (57). 
We adopt the latter approach. The prescriptions both assume that spin noise originates from a 
stationary process. This assumption may be appropriate if such spin noise originates from su- 
perfluid turbulence within the neutron star (/5) or fluctuations within the pulsar magnetosphere, 
but it may be a poor approximation if such noise arises from e.g. switches between metastable 
equilibria of the magnetosphere (/6, 58, 59). Most proposed mechanisms predict identical spin 
noise in the radio and gamma-ray bands. 
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Effects of the Ionized Interstellar Medium The ionized interstellar medium (IISM) is tur- 
bulent (60) and contains structures that act as lenses (67). Radio waves from a pulsar encounter 
a continually changing electron density and bend proportionally to the inverse square of the 
observing radio frequency v~*. The wide range of observational consequences includes dis- 
persive delays, strong intensity variations (scintillation), broadening of pulses due to multipath 
scattering, and higher-order effects, such as apparent position shifts (79, 62, 63). 

The line of sight through the IISM changes constantly due to the relative motion of Earth and 
the source. The most noticeable effect is scintillation: frequency- and time-dependent variations 
in the received intensity. Scintillation can render some observations useless if the pulsar is too 
faint, and it shifts the effective observing frequency by enhancing or depressing the signal over 
the observational band, magnifying the effects of DM uncertainty discussed below. 

The dominant effect of the ITSM on residuals is the dispersive delay T(t) x DM(t)/v?, 
which changes measurably on timescales of days to weeks (64). Low-frequency measurements 
yield higher DM precision, and since their inception, PTAs have monitored DM(t) using mul- 
tiple receivers covering a wide range of frequencies (e.g. 350 MHz to 4GHz). Wider band- 
width receivers (65) and high-cadence observations by low-frequency observatories (/8, 66) 
also improve DM measurements. However, the apparent DM also depends on v because lower- 
frequency radio waves scatter through a larger volume of the IISM, so even precise DM(t) 
measurements cannot fully correct (t) at other radio frequencies (17, 67). 

DM variation contributes potentially hundreds of degrees of freedom to a pulsar timing 
model. One approach uses DMX parameters (/3) to tabulate DM measurements from TOAs 


within discrete time segments. Another makes use of a constrained stochastic model for DM(t) 
(see spin noise above), which has fewer effective degrees of freedom but can only model sta- 
tionary DM variations, which is generally insufficient to capture observed variations (61,64, 66). 
These different approaches to DM modeling can produce discrepant results: NANOGrav, using 
DMX (9), and the PPTA, using the stochastic model (/4), find marginally different values for 
spin noise for pulsars common to the two PTAs (see below). 

Other IISM effects can be approximated using corrections of the form v^, e.g. v~* for 


multipath scattering, and other powers for higher-order effects (69), though only v? corrections 
are in widespread use. Scattering can be corrected precisely for a few bright MSPs using cyclic 
spectroscopy (70). A more universal method to incorporate such corrections is using a pulse 
portraiture, a model of a pulse profile over typically a factor of two or more of bandwidth (71). 
This portrait attempts to track the intrinsic variation in the pulse shape with frequency, and then 
is convolved with models of the IISM, e.g. a v~ kernel to account for a changing DM, a v^! 
kernel for altered scattering, etc. 

In summary, the monitoring, measurement, and interpolation of IISM-induced residuals is 


an observationally and computationally expensive endeavor and a source of unmodeled error. 
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Radio PTAs have managed to reduce these effects to the S1 ps level. None of the IISM effects 
described above, or solar wind effects described below, affect gamma-ray timing, which is also 
immune to the eclipses and DM variations affecting compact interacting binaries (72). 


The solar wind The solar wind makes a small contribution (typically 1074-1073 pc em 3) to 
the total DM for any pulsar, but it varies daily due to solar activity (20, 73) and annually, as the 
apparent pulsar-Sun angular separation changes. Pulsars near the ecliptic, e.g. PSR J00304-0451, 
are the most strongly affected. Existing models and observational strategies are insufficient to 
capture the contribution of the solar wind to pulsar timing residuals (73), so most PTA data 
analyses either combine solar wind variations with DM variations or employ a static model. 
Because the coherent structure of the solar wind subtends large portions of the sky, uncorrected 
solar wind variations can contribute correlated noise to pulsars (74). 


Pulse profile variations Pulse profile variations induce apparent variations in TOAs and can- 
not be easily mitigated. Observed variations in PSR J1643— 1224 (69) and PSR J0437—4715 
(48) have been attributed to a transient magnetospheric reconfiguration, while alterations in the 
pulse profile of PSR J17134-0747 may be associated with rapid apparent DM variations and 
recovery (75). No pulse profile variations have been observed in gamma-ray MSP profiles, 
though correlated variations in intensity and pulse shape have been observed for the young 
PSR J2021--4026 (76). A stable gamma-ray pulse profile is expected due to the likely origin 
of pulsed gamma-ray emission is from synchro-curvature radiation precipitated by large-scale, 
stable electric fields in the outer magnetosphere or current sheet (77). 


Discontinuities Any observation of a pulsar is referenced to international time standards, a 
well-established process with typical errors expected to be «10 ns. However, instrumenta- 
tion changes can alter the signal propagation time through the observing system, shifting the 
measured pulse phase relative to its true value. These phase shifts must be measured from or 
modelled in the data, incurring potential systematic error or adding typically 10—20 degrees of 
freedom in timing models (48). Fitting for these shifts acts as a high-pass filter, reducing sen- 
sitivity to low-frequency signals. Although not a generic property of gamma-ray pulsar timing, 
the LAT has a stable GPS clock and the LAT itself has been operating almost continuously with 
nearly-constant instrumental properties (23). Any calibration errors are static, resulting in static 
biases but no time-domain noise. 


Implications for gamma-ray analysis In summary, mitigating the noise sources for radio 
observations requires multi-frequency data, large fractional bandwidths, and homogeneous and 
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regular monitoring, a substantial practical challenge. In contrast, gamma-ray data only re- 
quire spin noise and Poisson noise models. This eases computational requirements and reduces 
systematic uncertainty. Due to continual all-sky monitoring, when a new MSP is discovered, 
archival LAT data can provide a full pulse timing history. The data span for each pulsar is 
uniform, ensuring that each pulsar is sensitive to the same spectrum of gravitational waves and 
enabling simple computational approaches. 


Data Preparation 


Photon data Properties of each gamma ray were determined by reconstructing the particle 
interactions in the LAT detector into measured quantities of incident direction, energy, and 
arrival time (22), with timestamping precision <300 ns (23). We began with a set of 127 MSPs 
with gamma-ray counterparts in 4FGL-DR2, the second release (78) of the fourth Fermi-LAT 
gamma-ray source catalog (79). For each pulsar, we downloaded data in the form of FT1 files 
(event lists) and FT2 files (tabulations of spacecraft position) from the Fermi Science Support 
Center (80) and processed it using the Fermi Science Tools v2.0.0 (87). We selected all data 
between 2008 Aug 04 (Modified Julian Day (MJD) 54682) and 2021 Jan 28 (MJD 59242) 
with reconstructed energy between 0.1 GeV and 10 GeV, measured zenith angle <100°, and a 
reconstructed incidence direction placing the photon within 3^ of the pulsar position. Using 
the 4FGL-DR2 sources models, we assigned each photon a weight (as in Equation S1) using 
the gtsrcprob Science Tool and restricted attention to events with w; > 0.05. This selection is 
intended to retain the great majority of the pulsar signal while eliminating background photons 
that increase computational costs. We used the PINT software package (82) to evaluate timing 
models and assign spin phase to photons. All raw and processed data, pulsar timing solutions, 
and software developed for this work are available in our data release (34). 


Ephemerides Of the 127 MSPs with 4FGL-DR2 spectral models, only 114 had suitable initial 
timing solutions, produced using data from the Nangay Radio Telescope (NRT), Green Bank 
Observatory (GBT), Arecibo Observatory (AO), the Parkes telescope/Murriyang (PKS), Jodrell 
Bank Observatory (JBO), and the Giant Metre-wave Radio Telescope (GMRT) (see Table S2). 
The parameters of these timing solutions were generally well constrained, but in some cases 
(noted in Table S2), the LAT data provided a more precise measurement. We re-fitted the 
timing model parameters using a maximum likelihood method, generally finding consistent 
results with the radio values. These parameters were allowed to vary in the photon-by-photon 
GWB analysis, but not in the TOA-based analysis. v and ù, the spin frequency and spin-down 
rate, were allowed to vary in both analyses. 
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Table S2: Properties of the 35 MSP ephemerides used for GWB analyses. Precise coordi- 
nates are provided in the data release (34). The observatory column indicates the primary source 
of data (acronyms defined in the text). “LAT” entries jointly fit the indicated parameter family 
with the GWB signal, e.g. position and proper motion (“astrometry”) or “binary” parameters 
such as the orbital frequency (FBO) and its time derivatives (FB1, FB2, ...). When these indi- 
cated parameters are fixed instead of fitted, the relative limit on the GWB drops from 1.0 to the 
value indicated in the fourth column. The final column indicates pulsars which use fewer than 
five frequencies to represent stochastic GWB or other red noise signals. The second, alternative 
model for PSR J1959+-2048 is discussed in the text. 


Name Observatory LAT Parameters Change | Harmonics 
PSR J0030+0451 | NRT - - - 
PSR J0034—0534 | NRT - - - 
PSR J0101—6422 | LAT binary,astrometry 0.95 - 
PSR J01024-4839 | NRT - - - 
PSR J0312—0921 | LAT binary,astrometry 0.86 3 
PSR J0340--4130 | GBT - - - 
PSR J04184-6635 | LAT astrometry 0.95 - 
PSR J0533+6759 | LAT astrometry 0.95 - 
PSR J0613—0200 | NRT - - 2 
PSR J0614—3329 | NRT - - - 
PSR J07404-6620 | NRT - - - 
PSR J1124—3653 | LAT binary,astrometry,FB1 0.92 — 
PSR J1231—1411 | NRT - - - 
PSR J1513—2550 | NRT - - 2 
PSR J1514—4946 | LAT binary,astrometry 0.98 - 
PSR J1536—4948 | LAT binary,astrometry 0.94 — 
PSR J1543—5149 | PKS - - 1 
PSR J1614—2230 | NRT - - - 
PSR J1625—0021 | LAT - - - 
PSR J1630+3734 | JBO - - - 
PSR J1741+1351 | NRT - - 2 
PSR J1810+1744 | LAT binary,astrometry,FB 1-8 0.73 - 
PSR J1816+4510 | LAT binary 0.97 - 
PSR J1858—2216 | LAT binary,astrometry 0.79 1 
PSR J1902—5105 | LAT binary,astrometry 0.94 — 
PSR J1908+2105 | LAT binary,astrometry 0.89 - 
PSR J19394-2134 | NRT - - - 
PSR J1959+2048 | LAT position 0.88 - 
PSR J1959+2048 | LAT position, binary, FB 1-10 0.46 - 
PSR J2017+0603 | NRT - - - 
PSR J20344-3632 | LAT astrometry 0.87 2 
PSR J2043--1711 | NRT - - - 
PSR J22144-3000 | NRT - - - 
PSR J2241—5236 | PKS+LAT binary 0.97 - 
PSR J2256—1024 | GBT+NRT+LAT | binary 0.91 - 
PSR J2302--4442 | NRT - E = 
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Solar system ephemeris In all cases, we used the DE421 Solar System ephemeris (83), which 
is commonly used in pulsar timing analyses and facilitates comparison with previous results. 
The choice of ephemeris can alter residuals at the ~100 ns level (84, 85), though see (/0), but 
the LAT data are not yet sensitive to such small effects. 


Pulse profile templates We modeled the pulse profile, f(¢), using 1-8 Gaussian components 
wrapped to the interval 0 < @ < 1, which enforces periodicity. We also fitted these using 
maximum likelihood, with brighter pulsars having more components. We held these templates 
fixed in the GWB studies described below. 


Photon-by-photon GWB Analysis 


Our photon-by-photon GWB analysis is implemented with maximum likelihood techniques. 
Degrees of freedom for spin noise or the GWB are incorporated into the log likelihood, 


log £ = Y log (wif (Olt AB) - 0— wi)) ~ 0.55 CS A. —0.5detC. (82) 
i kl 


The additional timing model parameters ( are the coefficients of the Fourier transform of a 
potential noise signal in the data, such that 


b(t, A, B) = e(t; A) + V a dba COS (2nk d ) + Pox41 sin (2i ) , (83) 


tobs obs 


with v the pulsar spin frequency and t_ps the time span of the data. The assumption that the 
noise follows a Gaussian random process is incorporated into the likelihood with the diagonal 
covariance matrix Cj; = da P (fk, Agwb)/2 which constrains the amplitudes 8. P(f,) is the 


power spectrum evaluated at the frequency fx = k/tobs (57). For a GWB, theory predicts 
2 —T=13/3 

P(f) = ne L i yr ?. In the analysis below we restricted 1 « k < 5 because, 

as also found by radio PTAs (9), we have verified that higher frequencies do not contribute 

significantly to constraints on the GWB. For some faint pulsars—typically those with additional 


free timing model parameters—the data are insufficient to constrain all five harmonics, so we 


used fewer frequencies as indicated in Table S2. 


MSP sample selection Only some MSPs are suitable for a GWB analysis. For each of the 
114 pulsars with initial radio timing solutions, we first estimated the white noise level, i.e. the 
typical amplitude of random fluctuations from the timing model, by averaging the measured 
amplitudes of 5 harmonics (6; in Equation S2). In the absence of other noise, this value is a 
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Figure S1: Timing properties of the 114 MSPs in the parent sample. The y axis shows the 
sensitivity, as estimated by the limiting amplitude of a white noise process, while the x axis 
gives the source significance, or brightness, as estimated by its total log likelihood, log £. Gray 
dashed lines indicate sample divisions, with the 29 orange stars indicating pulsars used in both 
analyses and the 6 purple triangles, those MSPs suitable for photon-by-photon analysis only. 
Gray circles indicate remaining MSPs. Four of these—PSR J1311— 1340, PSR J1555—2908, 
PSR J2215+5135, and PSR J2339—0533—lie within the selection region but are excluded as 
discussed in the text. 
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proxy for sensitivity. We also calculated the total significance of the pulsed signal via log £, 
which is a proxy for the number of independent TOAs that can be obtained for each pulsar. 
These two values, shown for the full sample in Figure S1, are correlated but with substantial 
scatter. At a given intensity, a pulsar with a narrower pulse or faster spin frequency has better 
timing precision. 

Faint MSPs cannot constrain the GWB but would increase computational complexity, so we 
required a white noise level < 500 us? yr^!, corresponding to a typical one-year TOA precision 
of ~16 us. We further set log £ > 150 for the photon-by-photon analysis to reduce problems 
with underconstrained parameters (see below). In addition to these requirements, we eliminated 
four pulsars. PSR J1311—3430 and PSR J1555— 2908 have strong spin noise, far in excess of 
a possible GWB signal, with sufficient amplitude to cause numerical issues in the algorithms 
described below. PSR J2215+5135 and PSR J2339—0533 have strong orbital period variations 
that require more model components than are supported by PINT. This selection narrowed our 
sample from 114 pulsars to 35. 


GWB Analysis For each pulsar we analyzed a potential GWB signal by evaluating the like- 
lihood following Equation S2, as a function of the timing model parameters A (v, 7, plus any 
LAT-constrained parameters as indicated in Table S2); the Fourier coefficients 8 representing 
the T = 13/3 GWB noise process; and the GWB amplitude Ayyp. Remaining timing model 
parameters were held fixed at radio values; the timing models in our data release (34) provide 
the exact degrees of freedom for each pulsar. 

A and P are nuisance parameters. To isolate A;,4, we marginalized over them by expanding 
the log likelihood as a quadratic form and solving the resulting Gaussian integrals; the likelihood 
surface for 8 is generally Gaussian (86). We did this for each pulsar, scanning Agwb over the 
range 10 ??-10- 9, The marginal log likelihoods obtained in this way, as a function of Agwb, are 
shown in Figure S2. With a uniform prior, they are equivalent to the logarithm of the posterior 
probability distribution. 

We obtained single pulsar limits on Agwb, listed in Table S4, by integrating the posterior 
probability density function until accumulating 9596 of it. In no single pulsar is there a strong 
detection of a GWB-like timing noise process, and in many cases the posterior probability 
distribution plateaus at Awb = 0. Joint limits on a common mode are obtained by multiplying 
the posteriors and integrating the resulting distribution. (This is the “factorized likelihood" 
approach (9).) The models for PSR J20434-1711 and PSR J2256— 1024 are somewhat improved 
with Aswb > 0, and when included in the common mode limit, they substantially increased it. 
We discuss these two pulsars in more detail below. 

To gauge the effect of non-GWB degrees of freedom (e.g. proper motion, binary parame- 
ters), we also obtained limits with those parameters fixed to their maximum likelihood values. In 
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Figure S2: The single-pulsar log likelihoods as a function of A, produced by the photon- 
by-photon method. Each likelihood was marginalized analytically over the timing model pa- 
rameters and Fourier coefficients of the noise process. The pulsars plotted with colored lines 
(see legend) have an individual limit Agwb < 1.5 x 107! or a peak value at Aswb > 0 along 
with a limit Ag, < 3.0 x 10713. These are the pulsars with the greatest influence. The other 
pulsars are shown in gray dotted lines. 
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Figure S3: The single-pulsar log likelihoods with an additional, numerically marginalized 
intrinsic spin noise process. Other content is as in as Figure S2. 


all cases, this improved the resulting limit, typically by 5-10% (see Table S2), and up to 15-25% 
for a few binaries. In the cases of PSR J1858—2216 (improved by 21%) and PSR J2034+3632 
(improved by 13%), the changes can be attributed to improved numerical stability due to the 
reduced degrees of freedom. Marked differences also occur for pulsars with strong variations 
in the orbital period, PSR J1810+1744 and PSR J1959+2048, indicating possible degeneracy 
between these parameters and the GWB parameters. Because the time scales associated with 
these processes are widely separated, the improved precision could also result from fitting many 
fewer degrees of freedom, viz. the 8 and 10 orbital frequency derivatives. 

Although we found no strong evidence for spin noise in our sample (see below), we also 
considered the case of per-pulsar spin noise along with a GWB signal. For each pulsar, we 
introduced a general red noise (RN) process with amplitude Ary and spectral index Ign, with 
Aprn defined identically to Agw» while In is allowed to take values other than 13/3. The power 
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spectral density Prn(f) is also defined identically to P4, ( f), and we evaluated the likelihood 
as above but now constraining the Fourier coefficients 8 by the summed power spectral density 
Py (f) + Pan(f). To marginalize over these two new nuisance parameters, we scanned over 
a grid of Ary and 1 < [pn < 7 and tabulated the resulting marginalized log likelihood as 
a function of Awb, again using a uniform prior on both Agwb and Arn. The log likelihoods 
in this case are shown in Figure S3. The peaks that appeared in Figure S2 have reduced or 
vanished, indicating that they are likely caused either by statistical fluctuations (i.e. they are low 
significance) or by weak timing noise processes with TrN 7 13/3. The resulting constraints on 
Aswb (Table S4) are generally poorer, except in cases where the log likelihood with an intrinsic 
RN model peaked at A,yp > 0. The most conservative limit, reported in the main text, was 
obtained using this approach. 


TOA-based GWB analysis 


The likelihood based-method described above takes advantage of the time resolution of the 
LAT data and, by construction, avoids potential systematic errors from reducing the full photon 
data to TOAs. However, TOA-based methods are computationally efficient, well-tested, and 
commonly used by radio PTAs, so we have implemented a parallel TOA-based analysis for 
comparison. 


TOA estimation The total log likelihood in Equation S1 can be interpreted as a measure of 
pulsation significance because f(@) = 1 and log Ldata = 0 in the absence of pulsation. log Ldata 
grows as data are accumulated, and once it surpasses a threshold (roughly 20) it is generally 
possible to reliably measure a phase shift 0¢ relative to the template. We converted dd to a TOA 
by choosing a reference time near the mid-point of the integration to such that ó(t9) mod 1 = 0 
according to a timing model, then iteratively determining ôt such that (to + 0;) mod 1 = dé. 
The resulting TOA, to + ôt, is in the Universal Time Coordinated (UTC)(GPS) time system 
because all timestamps are referenced to the on-board GPS clock (39). The shape of Z (09) 
becomes more Gaussian as the peak value increases, so a higher log £ threshold reduces the 
systematic uncertainty in TOA estimation (25) at the expense of a longer integration period. 
To minimize non-Gaussian effects, we required log £ > 300 to include a pulsar. For the 29 
suitable MSPs, we used the timing models resulting from the photon-by-photon analysis (with 
any LAT-optimized parameters) and computed TOAs with a cadence of 2, 1.5, and 1 TOAs per 
year, yielding a total of 25, 19, and 12 TOAs per pulsar, and analyzed (below) each of the data 
sets before choosing a single representative cadence for the final GWB analysis. 

We first applied TEMPONEST (26), which can constrain a variety of single-pulsar noise 
models and produce single-pulsar GWB limits via nested sampling methods. In particular, we 
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Table S3: Prior ranges for the noise models used in TEMPONEST and ENTERPRISE. The 
amplitudes are scaled such that the resulting power spectral density follows Equation 2. 


Parameter Prior ranges 
EFAC 0.1 to 5 
logio (EQUADs !) | —9 to —5 
logio ARN —18 to —10 
TRN 0 to 7 

logio Agwb —18 to —9 


tested for the presence of intrinsic red noise (RN) and excess white noise (WN). RN models, 
like the GWB, were assumed to be stationary processes with a power-law power spectral den- 
sity. WN models were implemented through the parameters EFAC and EQUAD (see above). In 
general, we expect no excess WN in the gamma-ray data. However, information is lost when 
reducing the full photon data to a single TOA and (assumed Gaussian) uncertainty. This could 
take the form of increased scatter of TOAs and appear as a WN process. 

Table S3 reports the priors on the noise model parameters. A slightly different range of 
spectral indices (0-7 vs. 1-7) is adopted between the photon-by-photon and TOA-based analy- 
ses, but this has little impact for inference on a steep spectrum process like the GWB. We used 
k = 12 frequency components fy to model low-frequency noise processes. Noise process am- 
plitudes follow linear-exponent priors, i.e. the prior probability density p(A) x 104, which is 
equivalent to sampling uniformly from log; A. All timing model parameters aside from v and 
v were held fixed. In all cases, we computed the Bayesian evidence for models with and without 
WN and RN, and the results are presented in Table S4. We observed that the preferred noise 
model for each pulsar excludes RN and WN processes for all pulsars except PSR J1959--2048 
and PSR J2241—5236. Both of these pulsars exhibit orbital period variations, and the excess 
white noise (detected with very modest Bayes factor ~7), could indicate unmodeled period 
variations. 

In general, we found that the results, including the inferred limit on the GWB, were consis- 
tent over the three cadences. For the faintest pulsars, compared to longer integrations, we expect 
the 2 yr! cadence to exhibit more systematic errors associated with the possibly-poor Gaussian 
approximation. This was the case for PSR J0533--6759, PSR J07404-6620, PSR J19394-2134, 
and PSR J2034--3632, for which models with additional WN were modestly preferred. With a 
1.5 yr! cadence, the preferred model for these pulsars required no additional WN. On the other 
hand, the highest cadence provides more information about high-frequency noise, so when pos- 
sible we adopted the 2 yr ^! cadence. The cadence chosen for each pulsar is indicated in Table 
84. 
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Using these preferred noise models and cadences, we next used ENTERPRISE (27) to search 
for GWB signals from each pulsar individually and with a correlation analysis. The priors on 
the noise parameters were the same (Table S3). Figure 2 shows that the TEMPONEST and 
ENTERPRISE limits for individual pulsars are consistent. In the multiple-pulsar analyses, we 
focused on the simpler common mode process, for which the analysis produced a joint pos- 
terior probability distribution for a fixed-amplitude, identical-spectrum noise process with an 
independent realization in each pulsar data set. We also performed a correlation analysis that 
assumed quadrupolar Hellings-Downs cross-correlation amplitudes (8). These results (see Ta- 
ble S5) were essentially identical to the common mode analysis. 


Results 


Comparison of methods The two codes ENTERPRISE and TEMPONEST provided single pul- 
sar limits which were consistent with each other, so we compared TOA-based and photon-by- 
photon approaches. Aside from computational aspects—the TOA-based methods are sampled, 
while the photon-by-photon method is analytic—there are two primary differences: the photon- 
based approach avoids the assumption of Gaussianity on TOA uncertainties, and it retains sen- 
sitivity to all timescales. We expect the photon-based approach to generally be more precise due 
to these advantages. The agreement between these two methods provides us with confidence in 
our photon-by-photon approach. The few exceptions (discussed below) stem from the funda- 
mental differences in the methods, and they have little impact on the final GWB limits because 
none of these pulsars contribute strongly to the sensitivity of the timing array. 

PSR J0340+4130 was the largest outlier, with a photon-by-photon limit nearly twice that 
of the TOA-based methods. The degrees of freedom in the two methods are identical, but 
the photon-by-photon log likelihood peaks at Ass, > 0, suggesting a real fluctuation or weak 
timing noise. We speculate non-Gaussianity may coincidentally distort the TOAs for this pulsar 
in such a way as diminish this possible noise and thus reduce the GWB limit. 

For PSR J1536—4948, the photon-based method also delivered a higher limit, probably 
because the photon-by-photon model includes free binary and astrometry parameters. The 
compact binary PSR J1810--1744 has extensive orbital period variations, and these degrees 
of freedom are likewise inaccessible to the TOA-based method. When we fixed these binary 
parameters, the photon-by-photon dropped by ~30% and agreed with the TOA-based results. 

PSR J1858—2215 and PSR J20344-3632 are both faint pulsars, so our models used a re- 
stricted number of harmonics for numerical stability (see above). This degraded the sensitivity 
of the photon-based analysis. 
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Table S4: Single pulsar limits on A, for 35 pulsars in our sample. These results use TEM- 
PONEST (TN in column 4), ENTERPRISE (ENT. in column 5) and the photon-by-photon method 
(columns 6 and 7). Pulsars with only photon-by-photon limits are indicated with an asterisk. 
Data for PSR J1959+2048 and PSR J2241—5236 favor a model with white noise, while all 
others favor no additional noise. Most pulsars can be analyzed with a 2 yr ? (182 day) cadence, 
while six pulsars require longer integrations (1.5 yr 1, 243 day) to produce reliable TOAs. The 
single-pulsar Awp limits are all 95% credible levels. Column 6 (photon-based limits) includes 
only the GWB in the noise model, while the limits in column 7 also include an intrinsic noise 
model (RN) for each pulsar, and the parameters of this noise model are marginalized. 


Pulsar Cadence | Noise model | TN ENT. Photon | Photon+RN 

TOA/yr (favored) | x10714 | x10714 | x10714 x10-14 
PSR J0030+0451 2 None 7.54 7.77 7.61 8.44 
PSR J0034—0534 2 None 13.40 13.39 15.56 18.00 
PSR J0101—6422 2 None 18.63 18.94 18.80 22.16 
PSR J0102+4839 2 None 39.29 38.90 38.90 38.76 
PSR J0312—0921* - - - - 21.57 27.86 
PSR J0340+4130 2 None 26.13 26.54 48.26 58.21 
PSR J0418+6635* - - - - 32.59 36.41 
PSR J0533--6759 L5 None 21.66 22.14 21.83 26.23 
PSR J0613—0200 2 None 22.57 21.59 21.66 25.80 
PSR J0614—3329 2 None 4.15 4.20 3.94 4.36 
PSR J0740+6620 1.5 None 15.76 16.62 17.24 20.55 
PSR J1124—3653 2 None 15.84 15.58 14.34 16.77 
PSR J1231—1411 2 None 2.19 2.30 2.71 3.54 
PSR J1513—2550* - - - - 25.57 43.22 
PSR J1514—4946 2 None 38.89 38.54 41.85 34.90 
PSR J1536—4948 2 None 12.30 11.69 14.36 15.38 
PSR J1543—5149* - - - - 98.02 1356.47 
PSR J1614—2230 2 None 9.08 9.23 8.14 9.76 
PSR J1625—0021 1.5 None 30.70 30.74 28.01 31.07 
PSR J1630+3734 2 None 7.28 7.27] 7.91 9.08 
PSR J17414-1351* - - - - 84.24 120.30 
PSR J1810+1744 2 None 14.31 14.95 18.54 21.35 
PSR J1816+4510 2 None 34.91 35.61 39.09 41.31 
PSR J1858—2216 2 None 31.51 30.46 110.83 1417.63 
PSR J1902—5105 2 None 11.50 11.38 11.87 15.06 
PSR J1908--2105* - - - - 41.80 47.97 
PSR J1939+2134 1.5 None 12.86 13.04 10.24 12.99 
PSR J1959+2048 2 WN 8.25 8.12 6.15 7.84 
PSR J2017+0603 2 None 16.63 16.60 17.59 20.25 
PSR J2034+3632 1.5 None 21.62 22.21 38.82 74.75 
PSR J2043+1711 2 None 13.38 13.80 13.96 15.07 
PSR J2214+3000 2 None 29.67 30.44 28.74 38.60 
PSR J2241—5236 2 WN 13.80 13.66 14.05 16.39 
PSR J2256—1024 1.5 None 12.98 13.23 13.84 13.31 
PSR J2302+4442 2 None 11.74 11.67 12.26 14.82 
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Table S5: 95% credible upper limits on A,,,;,/10~'* from the combined samples. The pul- 
sars corresponding to each subset are specified in the text and are ranked by their single-pulsar 
GWB upper limits. The “Full” rows indicate the total sample for the two methods, 29 pulsars 
common to TOA-based and photon-by-photon, and 35 to photon-by-photon only. The final row 
includes the Hellings-Downs (HD) spatial correlations predicted for quadrupolar gravitational 
radiation (8). 


Subset ENTERPRISE | ENTERPRISE | Photon | Photon 
with RN with RN 
Best 2 1.89 1.94 1.84 2.02 
Best 3 1.71 1.74 1.50 1.66 
Best 9 1.15 1.19 1.06 1.16 
Full 29 1.12 1.04 1.14 1.06 
Full 35 — — 1.14 1.05 
Full 29 w/HD 1.06 — — — 


Combined limits and scaling We combined the single-pulsar limits discussed above to obtain 
an overall limit on the GWB amplitude. To assess the dependence of the total limit on partic- 
ular pulsars, we obtained results for both the full pulsar timing array and for some subsets. 
The best nine pulsars in ascending order of single pulsar upper limits are PSR J1231—1411, 
PSR J0614—3329, PSR J1959+2048, PSR J0030+0451, PSR J1630+3734, PSR J1614—2230, 
PSR J1939+2134, PSR J1902—5105 and PSR J2302+4442. These pulsars form the subsets 
listed in Table S5. For both methods, the limits steadily improved with additional pulsars. In 
the case of the TOA-based approach, the limit also improved with the inclusion of the spatial 
correlation information predicted for a GWB (8). The two full-array limits were nearly identi- 
cal, while the tightest constraint was Agwb < 1.0 x 10 14. 

As discussed above, the limit in the photon-by-photon case can be degraded by includ- 
ing pulsars whose A,wt posteriors peak at Agsy > 0, which can happen even in the absence 
of a GWB signal due to statistical fluctuations. For completeness, however, we considered 
limits computed when removing two MSPs with the strongest signals, PSR J2043+1711 and 
PSR J2256—1024 (which exhibits modest orbital period variations). This is justifiable if the 
log likelihood peak is due to some deficiency in the timing solution. When removing these two 
pulsars, we obtained limits of 9.8 x 10715 and 10.3 x 10^ P. 

We additionally considered any systematic effects from our choice of the DE421 Solar Sys- 
tem ephemeris, by performing the same limit calculation using a perturbative Bayesian model- 
ing software, BAYESEPHEM (87). This allowed us to model uncertainties in the Solar System 
planetary masses and orbital parameters while simultaneously constraining other pulsar noise 
parameters and the GWB. The ENTERPRISE results, obtained with and without BAYESEPHEM, 
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Table S6: 95% credible upper limits on 4,,4,/10 !* from the combined samples using 
BAYESEPHEM. Here, we compare results obtained with ENTERPRISE, in both cases using 
spatial correlations (8) but with and without BAYESEPHEM. 


Subset | ENTERPRISE | ENTERPRISE 
with BAYESEPHEM 


Best 2 1.90 1.81 
Best 3 1.70 1.64 
Best 9 1.17 1.16 
Full 29 1.06 1.08 


are reported in Table S6. The limits with BAYESEPHEM were slightly tighter for the few-pulsar 
subsets, which follows from the increased degeneracy between the GWB quadrupolar signature 
and dipolar-like effects from Solar System ephemeris errors (74). The full array results, on the 
other hand, contained enough information to separate the two signals. The agreement between 
our limits fixed at DE421 and those obtained with BAYESEPHEM indicates its use does not bias 
our results. 

In summary, we obtained an upper limit of Agwb < 10 ^, with two independent methods 
yielding both expected scalings and consistent values. The limit was unaffected by the Solar 
System ephemeris (DE421) we chose, and to particular realizations of noise modeling, e.g. the 
inclusion of per-pulsar spin noise. 


Other GWB models Using the photon-by-photon method, we calculated 95% upper limits on 
more general power-law GWB models, shown in Figure 3. For each value of a, we computed 
a photon-based limit as described above, without including any per-pulsar spin noise. The 
resulting scaling is approximately exponential, and the scaling largely reflects the fact that the 
limit is dominated by frequencies around 2/t,p;. Limits are slightly less constraining for models 
with more high-frequency power due to the limited number of Fourier coefficients used in the 
analysis. 


Spinnoise There was little evidence for intrinsic spin noise in any of the pulsars in our sample 
(Table S4 and Figure S2). PTAs have measured such intrinsic timing noise for some of them, but 
the amplitudes are generally below the single-pulsar sensitivity. PSR J19394-2134 has intrinsic 
timing noise evident over the decades-long data sets (45). We detect no spin noise and obtained 
limits somewhat below the values reported by radio PTAs (see below). 

Although we saw no strong evidence for spin noise in any pulsar, PSR J2043--1711 and 
PSR J2256— 1024 had posterior probability distributions which peaked at As > 0, but at suf- 
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Table S7: Published radio PTA results on the GWB and GWB-like signals. Dates given are 
fractional year of publication. Upper limits (u.l.) and central values are given in the *A,,4" 
column while the “Range” column gives central ranges of GWB-like candidate signals. These 
values are plotted in Figure 1. 


Label Year Reference | Aswb Range | Note 
x10715 | x10715 
PPTA 2006 2006.96 (88) 11 — | 95% u.l. 
PPTA 2013 2013.79 (89) 2.4 — | 95% u.l. 
PPTA 2015 2015.71 (90) 1.0 — | 95% u.l. 
PPTA 2021 2021.54 (10) 2.2 | 1.9-2.6 | 68% range 
EPTA 2011 2011.54 (55) 6.0 — | 95% u.l. 
EPTA 2015 2015.88 (91) 3.0 — | 95% u.l. 
EPTA DR2 2021.96 (11) 2.3-3.7 | 5-95% credible region 
NANOGrav 5-yr 2013.04 (92) 7.0 — | 95% u.l. 
NANOGrav 11-yr | 2018.38 (85) 1.45 — | 95% u.l. 
NANOGrav 12.5-yr | 2020.96 (9) 1.92 | 1.4-2.7 | 5-95% credible region 
IPTA DRI 2016.38 (45) 1.7 — | 95% u.l. 
IPTA DR2 2022.21 (12) — | 2.0-3.6 | 5-95% credible region 


ficiently low values of A,y that the combined photon-by-photon limit increases. The posterior 
peaks decreased when we included a spin noise model and marginalized over its parameters, in- 
dicating the possible noise is likely more consistent with a different spectral index. In general, 
the single-pulsar Agwb limits increased when including these extra degrees of freedom (Table 
$4), but decreased for pulsars with possible intrinsic RN. The total limit on A5, is slightly 
lower when including spin noise, but both values were within 10% of each other. 


Comparison to radio measurements Radio PTAs have published spin noise models for 
many MSPs, particularly two recent sets based on the NANOGrav 12.5-yr data set (/3) and the 
PPTA data release 2 (/4, 48) containing measurements for 14 and 10 MSPs, respectively. Both 
contain treatments for the time-varying DM (NANOGrav use DMX and the PPTA use a station- 
ary process model), so these should be unbiased measurements of the spin noise. Of these pul- 
sars, only 3 are in common with our y-ray MSP sample, PSR J0030+0451, PSR J0613—0200, 
and PSR J1939+2134. For these pulsars, the radio estimates for spin noise do differ, with the 
PPTA analysis finding steeper spectra for the two common pulsars in our sample. 

We computed 95% upper limits on a spin-noise process for each pulsar using the radio- 
PTA measured spectral index I', with results listed in Table S8. Aside from the changed value 
of I°, the procedure is identical to the photon-by-photon method we used to compute single- 
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Table S8: Comparison between measured radio PTA spin noise amplitudes and 95% Fermi 
PTA upper limits. The spin noise power spectra are given in s? yr ! evaluated at f = 1/yr. 
Fermi limits are computed similarly to GWB amplitude limits except using the value of I° 
measured in each row rather than [ = 13/3. The last columns give the ratio of the Fermi upper 
limit to the radio PTA power. 


Pulsar PTA D P(f) PTA P(f) Fermi Ratio 
PSR J0030+0451 | NANOGrav | 63 | 9.0 x 10718 | 8.3 x 10-5 | 0.9 
PSR J0613—0200 | NANOGrav | 2.1 | 1.5 x 10714 | 2.0 x 10713 13 
PSR J0613—0200 PPTA | 4.2 | 2.5 x 10716 | 5.4 x 10715 21 
PSR J1939+2134 | NANOGrav | 3.3 | 9.8 x 10719 | 6.6 x 107! 0.7 
PSR J1939+2134 PPTA | 5.4 | 1.8 x 10716 | 1.1 x 107% | 0.6 


pulsar GWB limits. The Fermi upper limits for PSR J0030+0451 and PSR J1939+2134 were 
lower than the measured radio values and statistically incompatible. This could be evidence 
for uncorrected IISM/solar wind effects leaking into the spin noise estimation in the radio data. 
As discussed above, radio emission from PSR J00304-0451 is particularly affected by the solar 
wind because its position passes near to the Sun each year. 

Although the number of MSPs here is too small to draw firm conclusions, there is some 
evidence for contamination of spin noise models (and, ultimately, GWB-like signals) by uncor- 
rected IISM effects. The most sensitive gamma-ray MSPs, PSR J1231— 1411, PSR J0614—3329, 
and others are not among the pulsars routinely monitored by radio PTAs, so we cannot compare 
further. 
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